use "${clean}newspaper_analysis.dta", clear
keep if north == 1 | south == 1

// recoding Alexandria, VA to the South (technically part of DC until 1844)
replace south = 1 if city == "Alexandria"
replace north = 0 if city == "Alexandria"
replace state = "Virginia" if city == "Alexandria"

// dropping DC as a border area
drop if state == "District Of Columbia"

// standardizing variables
egen county_surban = std(county_urban_percent)
egen county_srugged = std(county_ruggedness)
egen county_spostoffice = std(county_postoffice_area)

// keeping the sample consistent
cap drop sample
gen sample = 1
foreach var of varlist $county_npcontrols {
	replace sample = . if `var' == .
}

// creating some variables
local cutyear = 1860

cap drop period
gen period = 0
replace period = 1 if year > `cutyear'
replace year = year - `cutyear'


// matrices to store the results for graphing
matrix define R = J(9, 3, 0)



* regressions
* ------------

xtset pub_id

// random effects model
xtreg singular c.year##i.period##i.north ${county_npcontrols} if sample == 1, re cluster(pub_id)

lincom _cons + 1.north // northern pre-war intercept
matrix R[1,1] = `r(estimate)'
matrix R[1,2] = `r(estimate)' + 1.96*`r(se)'
matrix R[1,3] = `r(estimate)' - 1.96*`r(se)'

local beta_re_n_i = `r(estimate)'
local beta_re_n_i : display %4.3f `beta_re_n_i'
local se_re_n_i = `r(se)'
local se_re_n_i : display %4.3f `se_re_n_i'	



lincom _cons // southern pre-war intercept
matrix R[2,1] = `r(estimate)'
matrix R[2,2] = `r(estimate)' + 1.96*`r(se)'
matrix R[2,3] = `r(estimate)' - 1.96*`r(se)'

local beta_re_s_i = `r(estimate)'
local beta_re_s_i : display %4.3f `beta_re_s_i'
local se_re_s_i = `r(se)'
local se_re_s_i : display %4.3f `se_re_s_i'	



lincom year + 1.north#c.year // northern pre-war slope
matrix R[5,1] = `r(estimate)'
matrix R[5,2] = `r(estimate)' + 1.96*`r(se)'
matrix R[5,3] = `r(estimate)' - 1.96*`r(se)'

local beta_re_n_s = `r(estimate)'
local beta_re_n_s : display %5.4f `beta_re_n_s'
local se_re_n_s = `r(se)'
local se_re_n_s : display %5.4f `se_re_n_s'	



lincom year // southern pre-war slope
matrix R[6,1] = `r(estimate)'
matrix R[6,2] = `r(estimate)' + 1.96*`r(se)'
matrix R[6,3] = `r(estimate)' - 1.96*`r(se)'

local beta_re_s_s = `r(estimate)'
local beta_re_s_s : display %5.4f `beta_re_s_s'
local se_re_s_s = `r(se)'
local se_re_s_s : display %5.4f `se_re_s_s'	



// fixed effects model
xtreg singular c.year##i.period##i.north ${county_npcontrols} if sample == 1, fe cluster(pub_id)

lincom year + 1.north#c.year // northern pre-war slope
matrix R[8,1] = `r(estimate)'
matrix R[8,2] = `r(estimate)' + 1.96*`r(se)'
matrix R[8,3] = `r(estimate)' - 1.96*`r(se)'

local beta_fe_n_s = `r(estimate)'
local beta_fe_n_s : display %5.4f `beta_fe_n_s'
local se_fe_n_s = `r(se)'
local se_fe_n_s : display %5.4f `se_fe_n_s'	



lincom year // southern pre-war slope
matrix R[9,1] = `r(estimate)'
matrix R[9,2] = `r(estimate)' + 1.96*`r(se)'
matrix R[9,3] = `r(estimate)' - 1.96*`r(se)'

local beta_fe_s_s = `r(estimate)'
local beta_fe_s_s : display %5.4f `beta_fe_s_s'
local se_fe_s_s = `r(se)'
local se_fe_s_s : display %5.4f `se_fe_s_s'	


* graphing 
* --------
svmat R

// temp variable to index estimates
cap drop temp*
gen tempn = _n if R1 != .

gen tempc = tempn > 4
label define tempclab 1 "Slopes" 0 "Levels", modify
label values tempc tempclab

gen tempnorth = 0
replace tempnorth = 1 in 1
replace tempnorth = 1 in 5
replace tempnorth = 1 in 8

gen tempfemodel = 0
replace tempfemodel = 1 in 8/9

drop if R1 == 0
 
twoway ///
	(scatter R1 tempn if tempnorth == 1 & tempc == 0,	msize(large) mlwidth(none) mcolor("${blue}") 	msym(D) yaxis(1 2)) ///
	(scatter R1 tempn if tempnorth == 0 & tempc == 0,	msize(large) mlwidth(none) mcolor("${red}")	msym(S) ) ///
 	(rcap R2 R3 tempn if tempnorth == 1 & tempc == 0 & tempfemodel == 0,   lwidth(thick) lcolor("${blue}") 	lpattern(dash) ) ///
	(rcap R2 R3 tempn if tempnorth == 0 & tempc == 0 & tempfemodel == 0,   lwidth(thick) lcolor("${red}") 	lpattern(dash) ) ///
 	, ///
	yline(0, lwidth(thick) lpattern(dash) lcolor(gs2)) ///
	text(.215 1 "{&beta}=`beta_re_n_i'", size(medlarge)) ///
	text(.185 1 "se=`se_re_n_i'", size(medlarge)) ///
	text(.215 2 "{&beta}=`beta_re_s_i'", size(medlarge)) ///
	text(.185 2 "se=`se_re_s_i'", size(medlarge)) ///
	legend(off) ///
	xlab(1 "North" 2 "South", labsize(large)) ///
	xtitle("Baseline", size(large) margin(t=3)) 	///
	ylab(none, axis(2) labsize(large)) ///
	ylab(0(.1).45, axis(1) grid glcolor(gs7%50)  labsize(large)) ///
	ytitle("") ///
	xscale(range(0.5 2.5)) ///
	subtitle("Singular Level in 1860", box bexpand bcolor(gs2) color(white) size(medsmall)) ///
	fxsize(100) fysize(100) ///
	name(intercepts, replace) nodraw
	
 
twoway ///
	(scatter R1 tempn if tempnorth == 1 & tempc == 1,	msize(large) mlwidth(none) mcolor("${blue}") 	msym(D) yaxis(1 2)) ///
	(scatter R1 tempn if tempnorth == 0 & tempc == 1,	msize(large) mlwidth(none) mcolor("${red}") 	msym(S)) ///
	(rcap R2 R3 tempn if tempnorth == 1 & tempc == 1 & tempfemodel == 0,   lwidth(thick) lcolor("${blue}") 	lpattern(dash)) ///
	(rcap R2 R3 tempn if tempnorth == 0 & tempc == 1 & tempfemodel == 0,   lwidth(thick) lcolor("${red}") 	lpattern(dash)) ///
	(rcap R2 R3 tempn if tempnorth == 1 & tempc == 1 & tempfemodel == 1,   lwidth(thick) lcolor("${blue}") 	lpattern(solid)) ///
	(rcap R2 R3 tempn if tempnorth == 0 & tempc == 1 & tempfemodel == 1,   lwidth(thick) lcolor("${red}") 	lpattern(solid)) ///
 	, ///
	yline(0, lwidth(thick) lpattern(dash) lcolor(gs2)) ///
	text(.00325 4.8 "{&beta}=`beta_re_n_s'", size(medlarge)) ///
	text(.0025 4.8 "se=`se_re_n_s'", size(medlarge)) ///
	text(.00325 6.2 "{&beta}=`beta_re_s_s'", size(medlarge)) ///
	text(.0025 6.2 "se=`se_re_s_s'", size(medlarge)) ///
	text(.0085 7.8 "{&beta}=`beta_fe_n_s'", size(medlarge)) ///
	text(.00775 7.8 "se=`se_fe_n_s'", size(medlarge)) ///
	text(.0085 9.2 "{&beta}=`beta_fe_s_s'", size(medlarge)) ///
	text(.00775 9.2 "se=`se_fe_s_s'", size(medlarge)) ///
	legend(off) ///
	xlab(5 "North" 6 "South" 8 "North" 9 "South", labsize(large)) ///
	xtitle("     Baseline                          Fixed Effects", size(large) margin(t=3)) 	///
	ylab(none, axis(1)) ///
	ylab(0(.002).01, axis(2) grid glcolor(gs7%50) labsize(large)) ///
	ytitle("", axis(2)) ///
	xscale(range(4 10)) ///
	subtitle("Slope of Time Trend in Singular Adoption: 1800-1860", box bexpand bcolor(gs2) color(white) size(medsmall)) ///
	fxsize(200) fysize(100) ///
	name(slopes, replace)  nodraw
 
graph combine intercepts slopes, scale(1.2) xsize(100) ysize(50)
graph export "${output}Fig2_np_prewar_ns_compare.pdf", as(pdf) replace


